{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 用scikit-learn做聚类分析\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x297eb125ba8>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAADxdJREFUeJzt3V+IXOd5x/HfsxvNKvK65CJb4VqW\nVWhZiAx1KuFmEJRZKw2uY9Jeymxy48Je1A4udWkrdFFyIXTl0Avnokss2pJtxFLHUFya1mw8NYbj\nOF5XTr1R1oRUdl2HCNMGeTHZ0e48vZhZj/7M7pxdzcx7njPfDyw7c/Tu+uFF/unlOe97xtxdAIA4\nxlIXAADYHYIbAIIhuAEgGIIbAIIhuAEgGIIbAIL5RJ5BZnZZ0oeSNiVtuPvxQRYFANheruBum3H3\nDwZWCQAgF1olABCM5Tk5aWb/Jen/JLmkv3H3+S5j5iTNSdL+/fuPHT58uM+lxtRsNjU2xr+PzEMH\nc9HBXHS8/fbbH7j7VJ6xeYP719z9fTP7VUkvSvqqu7+83fjp6WlfXV3NXXCZ1et11Wq11GUkxzx0\nMBcdzEWHmS3nvX+Y6586d3+//f2KpOclPbD38gAAt6NncJvZHWZ259ZrSV+Q9NagCwMAdJdnV8lB\nSc+b2db4f3D37w60KgDAtnoGt7v/VNJvDaEWAEAO3M4FgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAI\nhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGEEaWZTp37pyy\nLEtdSlJ5ProMAJLLskwnT55Uo9FQpVLR0tKSqtVq6rKSYMUNIIR6va5Go6HNzU01Gg3V6/XUJSVD\ncAMIoVarqVKpaHx8XJVKRbVaLXVJydAqAZBElmVaWFjQxMTEji2PLMtUr9dVq9W0tLT08etRbZNI\nBDeABLb61evr61pYWNi2X92tr3369OkEFRcLrRIAQ7fVr242mzv2q+lrd0dwAxi6rX712NjYjv1q\n+trd0SoBMHTValVLS0s6f/68HnvssW371Vvj6GvfiOAGkES1WtX6+nrPMK5WqwT2TWiVAEAwBDcA\nBENwA0AwBDcABENwA0AwBDeAQuNRrrdiOyCAwuJRrt2x4gZQWBx5747gBjAQ/WhxcOS9u9ytEjMb\nl/S6pP9x90cGVxKA6PrV4uDIe3e7WXE/KenSoAoBUB65WhwLC/rcqVPS2Jh05Ii0sND1d1WrVZ0+\nfZrQvk6u4DazQ5K+KOmbgy0HQBn0bHEsLEhzc9r/859L7tI770hzc9uGN25k7t57kNk/Sjon6U5J\nf9atVWJmc5LmJGlqaurY4uJin0uNaW1tTZOTk6nLSI556BiVuVhZWdHFixd1//336+jRozf82edO\nnWqF9k1+efCgXr1wYVglFsrMzMyyux/PM7ZncJvZI5Iedvc/NrOatgnu601PT/vq6mreekttqzc3\n6piHDuZCrfZIt+wxk5rN4ddTAGaWO7jztEpOSPqSmV2WdEHSg2b2rduoD8CoO3x4d9dxg57B7e6n\n3f2Qux+RdErS99z9ywOvDEB5nT0rHThw47UDB1rX0RP7uAEM3+ysND+vXx482GqP3HuvND/fuo6e\ndnXk3d3rkuoDqQTAaJmd1at3302/fw9YcQMjjAc4xcRDpoARxQOc4mLFDYyofjzAiRV7Gqy4gRG1\ndbpxa8W9214zK/Z0CG5gRN3uA5y6rdgJ7uEguIERVq1W9xy2t7tiTynLstBPHCS4AexJ1EeulqHF\nQ3AD2LPbWbGnUoYWD7tKAIyUMnyqDituACMlaovnegQ3gJETscVzPVolAEplFA4FseIGUBpl2DGS\nBytuAJLKsVLtxzH+CFhxAyjNSjXyoaDdILgBlGJvs1SOHSN5ENwASrVSjb5jJA+CG8DIrFTLguAG\nIGk0Vqplwa4SAAiG4AaAYAhuAAiG4AaAYAhuAAiG4AaAYAhuAAiG4AaAYAhuAAiG4AaAYAhuAAiG\n4AaAYAhuAAiG4AaAYAhuAAimZ3Cb2X4ze83M3jSzFTP72jAKAwB0l+eDFNYlPejua2a2T9IrZvYv\n7v7qgGsDAHTRM7jd3SWttd/ua3/5IIsCAGzPWrncY5DZuKRlSb8h6Rvu/hddxsxJmpOkqampY4uL\ni30uNaa1tTVNTk6mLiM55qGDuehgLjpmZmaW3f14nrG5gvvjwWafkvS8pK+6+1vbjZuenvbV1dXc\nv7fMtj58ddQxDx3MRQdz0WFmuYN7V7tK3P0XkuqSHtpDXQCAPsizq2SqvdKWmX1S0ucl/XjQhQEA\nusuzq+QuSX/X7nOPSVp09xcGWxYAYDt5dpX8UNJnh1ALACAHTk4CQDAENwAEQ3ADQDAENwAEQ3AD\nQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAE\nNwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAE\nQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDA9g9vM7jGzl8zskpmtmNmTwygMANDdJ3KM2ZD0lLu/YWZ3\nSlo2sxfd/UcDrg0A0EXPFbe7/8zd32i//lDSJUl3D7owAEB35u75B5sdkfSypPvc/epNfzYnaU6S\npqamji0uLvavysDW1tY0OTmZuozkmIcO5qKDueiYmZlZdvfjecbmDm4zm5T075LOuvt3dho7PT3t\nq6uruX5v2dXrddVqtdRlJMc8dDAXHcxFh5nlDu5cu0rMbJ+k5yQt9AptAMBg5dlVYpKelXTJ3b8+\n+JIAADvJs+I+Iekrkh40s4vtr4cHXBcAYBs9twO6+yuSbAi1AABy4OQkAARDcANAMAQ3AARDcANA\nMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcBdMlmU6\nd+6csixLXQqAgur5QQoYnizLdPLkSTUaDVUqFS0tLalaraYuC0DBsOIukHq9rkajoc3NTTUaDdXr\n9dQlASgggrtAarWaKpWKxsfHValUVKvVUpcEoIBolRRItVrV0tKS6vW6arUabRIAXRHcBVOtVgls\nADuiVQIAwRDcbWzDAxAFrRKxDQ9ALKy4xTY8ALEQ3GIbHoBYaJWIbXgAYiG429iGByAKWiUAEAzB\nDQDBENwAEAzBDQDBENx9xOlLAMPArpI+4fQlgGHpueI2s/NmdsXM3hpGQVFx+hLAsORplfytpIcG\nXEd4nL4EMCw9WyXu/rKZHRl8KbEN6vRllmWc6ARwA3P33oNawf2Cu9+3w5g5SXOSNDU1dWxxcbFP\nJca2tramycnJPf3sysqKnnrqKV27dk379u3T008/raNHj/a5wuG4nXkoG+aig7nomJmZWXb343nG\n9u3mpLvPS5qXpOnpaadV0LK1Wt7RwoJ05oz07rvS4cPS2bPS7KyyLNPGxoaazaY2NjZ09erVsC2Y\nXPMwIpiLDuZib9hVktrCgjQ3J330Uev9O++03qvTN9/aqcJfcAASwZ3emTOd0N7y0UfSmTOqXr7M\nUwsB3KJncJvZtyXVJH3azN6T9Ffu/uygCxsZ77674/W9PrWQm5pAeeXZVfLoMAoZWYcPt9oj3a7v\n0fz8vJ544gltbm5qYmKCw0BAyXDkPbWzZ6UDB268duBA6/oeZFmmxx9/XNeuXVOz2dT6+jqHgYCS\nIbhTm52V5uele++VzFrf5+db1/egXq+r2Wx+/H58fJybmkDJcHOyCGZn9xzUN6vVapqYmND6+rrG\nxsb0zDPP0CYBSobgLhk+PxMoP4K7rSi7MPpRB5+fCZQbwa3iPJK1KHUAKDZuTqo4j2QtSh0Aio3g\nVnEeybrbOvjEHWA00SpRcW7o7aYO2irA6CK424pyQy9vHd3aKkWoH8Dg0SoJqijtHQDDx4o7qLxt\nlaJscwTQPwR3YL3aKvTBgXKiVVJibC8EyongLjH64EA50SopsaJscwTQXwR3yRVlmyOA/qFVAgDB\nENwAEAzBDQDBENwAEAzBDQDBENwAEAzBDQDBENwAEAzBDQDBENwAEAzBDQDBENwAEAzBDQDBENwA\nEAzBDQDBENwAEAzBDQDBENwAEEyu4Dazh8xs1cx+YmZ/OeiiAADb6xncZjYu6RuSfl/SZyQ9amaf\nGXRhAIDu8qy4H5D0E3f/qbs3JF2Q9AeDLQsAsJ08n/J+t6T/vu79e5J+5+ZBZjYnaa79dt3M3rr9\n8krh05I+SF1EATAPHcxFB3PRMZ13YJ7gti7X/JYL7vOS5iXJzF539+N5iygz5qKFeehgLjqYiw4z\nez3v2Dytkvck3XPd+0OS3t9tUQCA/sgT3D+Q9Jtm9utmVpF0StI/DbYsAMB2erZK3H3DzJ6Q9K+S\nxiWdd/eVHj8234/iSoK5aGEeOpiLDuaiI/dcmPst7WoAQIFxchIAgiG4ASCYvgY3R+NbzOy8mV1h\nL7tkZveY2UtmdsnMVszsydQ1pWJm+83sNTN7sz0XX0tdU2pmNm5m/2FmL6SuJSUzu2xm/2lmF/Ns\nC+xbj7t9NP5tSb+n1hbCH0h61N1/1Jf/QCBm9ruS1iT9vbvfl7qelMzsLkl3ufsbZnanpGVJfzii\nfy9M0h3uvmZm+yS9IulJd381cWnJmNmfSjou6Vfc/ZHU9aRiZpclHXf3XIeR+rni5mh8m7u/LOl/\nU9dRBO7+M3d/o/36Q0mX1DqNO3K8Za39dl/7a2R3B5jZIUlflPTN1LVE08/g7nY0fiT/B0V3ZnZE\n0mclfT9tJem0WwMXJV2R9KK7j+xcSPprSX8uqZm6kAJwSf9mZsvtx4fsqJ/BnetoPEaTmU1Kek7S\nn7j71dT1pOLum+5+v1onkB8ws5FspZnZI5KuuPty6loK4oS7/7ZaT2F9vN1u3VY/g5uj8eiq3c99\nTtKCu38ndT1F4O6/kFSX9FDiUlI5IelL7d7uBUkPmtm30paUjru/3/5+RdLzarWet9XP4OZoPG7R\nviH3rKRL7v711PWkZGZTZvap9utPSvq8pB+nrSoNdz/t7ofc/YhaWfE9d/9y4rKSMLM72jfuZWZ3\nSPqCpB13pPUtuN19Q9LW0fhLkhZzHI0vJTP7tqRM0rSZvWdmf5S6poROSPqKWiuqi+2vh1MXlchd\nkl4ysx+qtdB50d1HehscJEkHJb1iZm9Kek3SP7v7d3f6AY68A0AwnJwEgGAIbgAIhuAGgGAIbgAI\nhuAGgGAIbgAIhuAGgGD+H/AP/bWDoznWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x297e90cb240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 构造数据\n",
    "cluster1 = np.random.uniform(0.5, 1.5, (2, 10))\n",
    "cluster2 = np.random.uniform(3.5, 4.5, (2, 10))\n",
    "\n",
    "# 顺序连接两个矩阵，形成一个新矩阵,所以生成了一个2*20的矩阵，T做转置后变成20*2的矩阵,刚好是一堆(x,y)的坐标点\n",
    "x = np.hstack((cluster1, cluster2)).T\n",
    "plt.figure()\n",
    "plt.axis([0, 5, 0, 5])\n",
    "plt.grid(True)\n",
    "plt.plot(x[:, 0], x[:, 1], 'k.')\n",
    "# 使用k-means做聚类\n",
    "kmeans = KMeans(n_clusters=2)\n",
    "kmeans.fit(x)\n",
    "plt.plot(kmeans.cluster_centers_[: ,0], kmeans.cluster_centers_[:, 1], 'ro') # 红色为重心点位置"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 肘部法则（现实情况会比以上构造的数据复杂很多，因为直接看根本不知道该聚几个类。）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 随机生成4个群\n",
    "cluster1 = np.random.uniform(0.5, 1.5, (2, 10))\n",
    "cluster2 = np.random.uniform(1.5, 2.5, (2, 10))\n",
    "cluster3 = np.random.uniform(1.5, 3.5, (2, 10))\n",
    "cluster4 = np.random.uniform(3.5, 4.5, (2, 10))\n",
    "x_1 = np.hstack((cluster1, cluster2))\n",
    "x_2 = np.hstack((cluster3, cluster4))\n",
    "x = np.hstack((x_1, x_2)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "k = np.arange(1, 10)\n",
    "meandistortions = []\n",
    "for k_ in k:\n",
    "    kmeans = KMeans(n_clusters=k_)\n",
    "    kmeans.fit(x)\n",
    "    # 求kmeans的成本函数值\n",
    "    meandistortions.append(sum(np.min(cdist(x, kmeans.cluster_centers_, 'euclidean'), axis=1)) / x.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x297eb1ff470>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHOhJREFUeJzt3XuUXFWZ9/Hvk046zSUEJwkQCdiO\n4j1IoEEiM6SbDHIdwosOIndZrjAovuFVVoYwo+iwmIwsdcIoAytclCgD5A0MCwF5BzDNRRpMB7kE\nIxgRNYCkEyTcJE06z/vHrprutHXvU7XPqfp91qrVp6pOqp466X7Ormfvs7e5OyIi0lzGxQ5ARESS\np+QuItKElNxFRJqQkruISBNSchcRaUJK7iIiTUjJXUSkCSm5i4g0ISV3EZEmND7WG0+dOtU7Oztj\nvb2ISCatXr16o7tPK7dftOTe2dlJf39/rLcXEckkM/ttJfupLCMi0oSU3EUkk/r6+li8eDF9fX2x\nQ0mlaGUZEZFa9fX1MXfuXAYHB2lvb+fee+9l9uzZscNKFbXcRSRzent7GRwcZGhoiMHBQXp7e2OH\nlDpK7iKSOuVKLt3d3bS3t9PW1kZ7ezvd3d2NDXAMGlVOUllGRFKlkpLL7Nmzuffee+nt7aW7u7vh\nJZm+vj6WLVsGwOmnn17x+zeynJRocjezNqAfeN7dj03ytUWkNRQquRRKgLNnz45SZ+/r66Onp4ct\nW7YAcO211xaNcbRKP1sSki7LLADWJvyaItJC0l5yySfovLfffrvimn8jP1tiLXczmwEcA1wCfCmp\n1xWR1hK75FJOPkHnW+4TJkyoOEk38rNZUgtkm9kKYDEwCTi/XFmmq6vLdYWqiGRRrTX3JJjZanfv\nKrdfIi13MzsW2ODuq82su8R+84H5AHvvvXcSby0i0nCx6v3VSKrmfghwnJk9B9wIHGZmPxy9k7sv\ndfcud++aNq3svDciIlKjRJK7uy9y9xnu3gmcBPzE3U9N4rVFRKR6uohJRKRBGjkfTuIXMbl7L9Cb\n9OuKiMTQ19eXyOiWRs+HoytURUSKSDIhN/ICJlBZRkSkqCQnKGv0xVlquYtIVEmVPeohn5DzLfex\nJORGX5yV2EVM1dJFTCKSL3ts2bKFcePGcfnllzN//vzYYW0nbSefhl7EJCJSi97eXrZs2cK2bdvY\ntm0b5557LjNnzhxzx2WSyTgLFywVouQuItF0d3czbtw4tm3bBsDQ0NCYOhq1QtMwdahKS9L6m+kw\ne/ZsLr/8ciZMmMC4ceOYOHHimOraWqFpmFru0nLUukuX+fPnM3PmzERKKUl2gGadkrtULG0dS7Vq\n9HhjKS+punbapwtuJCV3qUgztXbVumtuWe0ATZqSu1SkmVq7at1JK1Byl4o0W2tXrTtpdkruUhG1\nduuvWfo0JB2U3KViau3WT5r7NHTSySYld5EK1DvBpbVPI80nHSlNyV2kjEYkuLT2aaT1pCPlKbmL\nlNGIBJfWPo20nnSkPCV3kTIaleDS2KeR1pOOlKcpf0UqoE5FSQtN+SuSoDS2qkVK0ayQIiJNSMld\nRKQJKbmLiDQhJXcRkSak5C4i0oSU3EVEmpCSu4hIE1JyFxFpQkruIiJNSMldRKQJKbmLiDQhJXeR\nMvr6+li8eDF9fX2xQxGpmCYOEykh7SsRabZKKUbJXaSENK9ElPYTj8SlsoxICfmFOtra2lK3ElGh\nE49InlruUpVWKwOkeSUiLYEnpSi5S8XqXQZI24ljZDyLFi2KHc6fSfOJR+JTcpeK1bP+nLb6cdri\nKUYrREkxqrlLxYrVn5MYKpi2+nHa4hGpllruUrFCZYCkWrhpqx+nLR6RaiWS3M1sL2AZsAewDVjq\n7pcl8dqSLqPLAMuWLeOtt97C3cdUqklb/Tht8YhUy9x97C9iNh2Y7u6PmtkkYDVwvLv/oti/6erq\n8v7+/jG/t8TT19dHT08PW7ZsAaC9vb1u48DT1tkqEouZrXb3rnL7JdJyd/cXgRdz26+Z2VpgT6Bo\ncpfs6+3tZevWrQCYGWeddVbdEnsWOjdF0iTxDlUz6wRmAY8UeG6+mfWbWf/AwEDSby0NNrKDtaOj\ng9NPP72m1ynXIavOTZHqJdqhamY7AzcD57n7q6Ofd/elwFIIZZkk31uSVUkZJIm6dCWtcnVuilQv\nseRuZhMIif16d78lqdeVxqumDFLLOOuRJ45Kxs6rc1OkekmNljHgGmCtu387ideUeBp5sdKSJUsq\napXrYh2R6iTVcj8EOA140sweyz12obvfmdDrSwPVswwy+sSxadOm1LfKNVJHsiip0TIPApbEa0l8\nSZVBCiXFQieONLfKNVJHskpXqEpBY024xZJi1urnaZ7PXaQUJXepi1JJMc0t9dE0UkeySsld6qKa\npJjmmnbWvmmI5CUy/UAtNP1A86skaaumLVKdhk4/0IrS3NpMi0rKL5XUtHWsRaqn5F4DtTaTU658\no2MtUhst1lEDzXWSnHxN++KLLy6YuHWsRWqjlnsNNIIiWaXKNzrWIrVRcq+BRlA0jo61SG00WkZE\nJEMqHS2jmruISBNSchcRaUJK7kWUWx1IRCTN1KFagMZWi0jWqeVegMZWi0jWKbkXMHLhZ42tFpEs\nUlmmAI2tFpGsU3IvIktzjouIjKayjIhIE1JyTzkNyRSRWqgs0yC1zEmuIZkiUisl9waoNUmXW8hC\ni1iISDFK7g1QyWpDhZSa7latehEpRcm9AWqdk7zUkMxaTxgi0hqaLrmnsVQxlnHzo4dk5j/flClT\ntIiFiBTVVMk9zaWKJMbNj/58S5YsYdOmTak6kYlIOmRuKGSpoYH1mhMmLcMRR3++TZs2sWjRIiV2\nEfkzmWq5l2uZ12O9zTR9Gyj3+dJYkhKRODKV3Mt1ItZjTpg0dVyW+nxpOgmJSHyZSu6VtMyTnhOm\nHt8GxqLY50vTSUhE4stUco8xW2NWZohM20lIROIyd4/yxl1dXd7f3x/lvZuVau4izc/MVrt7V7n9\nMtVyl9I0TbGI5EVruZvZAPDbGv/5VGBjguHUg2Icu7THB4oxKWmPMU3xvcvdp5XbKVpyHwsz66/k\na0lMinHs0h4fKMakpD3GtMdXSOYuYhIRkfKU3EVEmlBWk/vS2AFUQDGOXdrjA8WYlLTHmPb4/kwm\na+4iIlJaVlvuIiJSgpK7iEgTSnVyN7MjzexpM1tnZhcUeP5MMxsws8dyt881OL5rzWyDma0p8ryZ\n2b/n4n/CzPZvZHwVxthtZptHHMOvNji+vcxspZmtNbOnzGxBgX2iHscKY4x9HDvM7Gdm9nguxq8X\n2Geimd2UO46PmFlnyuKL+vc8Io42M/u5md1e4Llox7Bq7p7KG9AG/Br4S6AdeBz40Kh9zgS+GzHG\nQ4H9gTVFnj8a+DFgwMHAIymMsRu4PeIxnA7sn9ueBDxT4P856nGsMMbYx9GAnXPbE4BHgINH7fN5\n4Mrc9knATSmLL+rf84g4vgT8Z6H/z5jHsNpbmlvuBwHr3P1Zdx8EbgTmRY5pO+5+P/ByiV3mAcs8\neBjY1cymNya6oIIYo3L3F9390dz2a8BaYM9Ru0U9jhXGGFXu2Lyeuzshdxs9WmIecF1uewUw18ws\nRfFFZ2YzgGOAq4vsEu0YVivNyX1P4Pcj7q+n8B/UJ3Nf1VeY2V6NCa1ilX6G2Gbnvi7/2Mw+HCuI\n3FfcWYRW3UipOY4lYoTIxzFXTngM2ADc7e5Fj6O7bwU2A1NSFB/E/3teAiwEthV5PuoxrEaak3uh\ns+HoM/2PgE533xe4h+EzalpU8hlie5QwV8VHge8At8YIwsx2Bm4GznP3V0c/XeCfNPw4lokx+nF0\n9yF33w+YARxkZh8ZtUvU41hBfFH/ns3sWGCDu68utVuBx9L2Nw2kO7mvB0aeuWcAL4zcwd03ufuW\n3N2rgAMaFFulyn6G2Nz91fzXZXe/E5hgZlMbGYOZTSAkzevd/ZYCu0Q/juViTMNxHBHLK0AvcOSo\np/7nOJrZeGAyEUp2xeJLwd/zIcBxZvYcoQx8mJn9cNQ+qTiGlUhzcl8F7GNm7zazdkLnxW0jdxhV\ndz2OUAtNk9uA03OjPQ4GNrv7i7GDGsnM9sjXDM3sIMLvxKYGvr8B1wBr3f3bRXaLehwriTEFx3Ga\nme2a294B+Bvgl6N2uw04I7f9KeAnnusZTEN8sf+e3X2Ru89w905CvvmJu586ardox7BaqZ3P3d23\nmtm5wP8jjJy51t2fMrN/Bvrd/Tbgf5vZccBWwtnzzEbGaGY3EEZJTDWz9cBFhI4i3P1K4E7CSI91\nwJvAZxsZX4Uxfgo4x8y2An8CTmrwL+shwGnAk7l6LMCFwN4jYox9HCuJMfZxnA5cZ2ZthBPLcne/\nfdTfyzXAD8xsHeHv5aSUxRf177mYFB3Dqmj6ARGRJpTmsoyIiNRIyV1EpAkpuYuINKFoHapTp071\nzs7OWG8vIpJJq1ev3ugVrKEaLbl3dnbS399f8f6XXgoHHgg9PcOPrVwJq1bBwoV1CFBEJIXM7LeV\n7JeZssyBB8KJJ4aEDuHniSeGx0VEZHupHec+Wk8PLF8Of/d3Ybu3N9wf2ZIXEZEgMy13CIn8Ax+A\nFStg7lwldhGRYjKV3FeuhKefhj32CK32ZctiRyQikk6ZSe75Gvvy5fDwwzBpEpx1Ftz+Z2uliIhI\nZpL7qlXDNfZ3vQtuvRXc4YILYFuxmZdFRFpUZpL7woXb19h7euBb34KnnoJLLokXl4hIGmUmuRey\nYAGceipcdJHKMyIiI2U6uZvB0qUwaxaccgo880zsiERE0iHTyR1ghx3gllugvR2OPx5eHb34mYhI\nC8p8cofQwbp8eWi5n3GGOlhFRJoiuUPoYP3mN8MoGnWwikira5rkDupgFRHJa6rkrg5WEZGgqZI7\nqINVRASaMLmDOlhFRJoyuYM6WEWktTVtcgd1sIpI62rq5K4OVhFpVU2d3EEdrCLSmpo+uYM6WEWk\n9bREcgd1sIpIa2mZ5A7qYBWR1tFSyV0drCLSKloquYM6WEWkNbRccgd1sIpI82vJ5A7bd7D+y7/E\njkZEJFktm9xhuIP1q1+FO+6IHY2ISHJaOrmP7GA9+WR1sIpI82jp5A7qYBWR5lQ2uZvZtWa2wczW\nFHnezOzfzWydmT1hZvsnH2Z9qYNVRJpNJS337wNHlnj+KGCf3G0+cMXYw2o8dbCKSDMpm9zd/X7g\n5RK7zAOWefAwsKuZTU8qwEZSB6uINIskau57Ar8fcX997rHMUQeriDSLJJK7FXjMC+5oNt/M+s2s\nf2BgIIG3Tp46WEWkGSSR3NcDe424PwN4odCO7r7U3bvcvWvatGkJvHV9qINVRLIuieR+G3B6btTM\nwcBmd38xgdeNSh2sIpJl48vtYGY3AN3AVDNbD1wETABw9yuBO4GjgXXAm8Bn6xVsoy1YAKtXhw7W\nWbPgmGNiRyQiUhlzL1ger7uuri7v7++P8t7V+NOf4K/+Ctatg1Wr4H3vix2RiLQyM1vt7l3l9mv5\nK1TLUQeriGSRknsF1MEqIlmj5F4hdbCKSJYouVdBV7CKSFYouVdBV7CKSFYouVdJHawikgVK7jVQ\nB6uIpJ2Se43UwSoiaabkPgYLFoT6+1e+sn0H68qVcOml8eISEVFyHwOz0GofPx5OPDGUaVauDNsH\nHhg7OhFpZWXnlpHSjjwSli2DU06B2bND/f2WW0LZRkQkFrXcE/CZz8Bpp8HLL8Mrr8A118Dvfhc7\nKhFpZUruCVi5Eu68ExYuDEMlly+H978fLrxQQyVFJA4l9zHK19iXL4dvfCN0rE6aFGaSXLwY3vte\nuOIK2Lo1dqQi0kqU3Mdo1aqQ2PM19p4eWLECDj88PPfBD8LnPw8zZ8Ltt0OkGZZFpMVoPvc6c4fb\nbgslm2eegcMOg299C/bbL3ZkIpJFms89Jcxg3jxYswa+8x14/HHYf38480x4/vnY0YlIs1Jyb5AJ\nE+Dcc8OKTuefDzfcAPvsE2aYfO212NGJSLNRcm+wXXcNV68+/XRo0V98cUjyV12lTlcRSY6SeySd\nnaH1/vDDYUTN/PlhKoO77oodmYg0AyX3yD72MXjgAbj55rAY91FHwRFHwBNPxI5MRLJMyT0FzOCE\nE+AXv4B/+7cwhHLWLPjc5+DFF2NHJyJZpOSeIu3tcN558Otfh5/LloV6/Ne/Dm+8ETs6EckSJfcU\nesc7wlj4tWvh6KPha1+D970Pvvc9GBqKHZ2IZIGSe4q95z3h6tef/hT23hvOOgsOOADuuSd2ZCKS\ndkruGfDxj8NDD8FNN4WJyA4/PLTon3oqdmQiklZK7hlhFiYoW7s2LO/30EOw777w938PL70UOzoR\nSRsl94yZOBG+/OXQ6frFL4a549/7XrjkEnjzzdjRiUhaKLln1JQpsGRJGD55+OHwT/8U5pBftiys\nBiUira2i5G5mR5rZ02a2zswuKPD8mWY2YGaP5W6fSz5UKWSffcKyfvffD9OnwxlnhPVbzz47zDU/\nkhbuFmkdZZO7mbUBlwNHAR8CPmNmHyqw603uvl/udnXCcUoZf/3XYSqD66+HjRth6dKwvut114Xn\ntXC3SGuppOV+ELDO3Z9190HgRmBefcOSWowbByefHCYl+9d/hfHjw9TCBxwAxx8P3/++Fu4WaRWV\nJPc9gd+PuL8+99honzSzJ8xshZntlUh0UpOODviHf4DnnoOuLnj00TCEct48OOigsHDIHXfA5s2x\nIxWReqkkuVuBx0Yv3/QjoNPd9wXuAa4r+EJm882s38z6BwYGqotUqrZmTUjwF1wAkyeHVn1HB1x2\nGRx7LPzFX4Tk/+Uvh9Wi/vjH2BGLSFLKLrNnZrOBr7n7Ebn7iwDcfXGR/duAl919cqnXbZVl9mIZ\nuXB3T8/29w8+GB55BHp74b77oK8PtmwJY+k/+lGYMwe6u0Mdf8qU2J9EREaqdJm98RW81ipgHzN7\nN/A8cBJw8qg3m+7u+fkLjwPWVhmvJKzQwt3Ll4fHe3pC8u7uDs+99Rb87Gch0ff2hs7Yyy4Lz82c\nGfabMwcOPRSmTWv8ZxGR6lW0QLaZHQ0sAdqAa939EjP7Z6Df3W8zs8WEpL4VeBk4x91/Weo11XJP\nr8HBcBLIt+x/+tPhC6Q+/OHhlv2hh8Luu8eMVKT1VNpyryi514OSe3a8/Tb094dEf9998OCD8Prr\n4bkPfGC4ZT9nThhrLyL1o+QudbN1axiBky/jPPhgGI0D4aKqkcl+xoyYkYo0HyV3aZihIXjsseEy\nzgMPwCuvhOfe857hMs6cOXDjjeFCqpHj7VeuDGWghQtjRC+SLUruEs3QUFgDNl/Gue++4WGWu+8e\nEv+CBWEZwV/9KkyZMLLzV0SKU3KX1Ni2LYy5z5dx7rlnuIwD4Ura3XYLwy5L3aZOHd7edVdoa0su\nxksv1TcKyYYkh0KKjMm4cWHu+X33DdMUb9sGX/gCXHklzJ0bpkfYtGn4tnbt8PbWrYVf0ywsR1jq\nBFDo1tFR+PUOPLD4dQEiWaTkLg13332wYgV85StwxRXwj/9YuCTjHlr4IxP/yNvGjcPbL7wATz4Z\ntkstJr7jjsVPAJ/+dJii4dhj4cc/hv/4jzBdg3s4mYhkicoy0lClrpxNqub+1lvw8svbJ/9SJ4ZN\nm0KfQLE/hY6OcEKYOjVcxJXfLvXYhAnJfBZQyUi2p7KMpFK5K2eT0NEB73xnuFVqaAh+9KOwCPm8\neWGO/LPPDol640YYGAg/N26E3/wm3C818drkycWTf6H7kyeH8lUhKhlJLdRyF6G2bxSDg8PfEEYm\n/5HbI+8PDIQ5fAppaxvuMyj0TWBgAL77XfjkJ+HWW0N/xdFHw047qWTUajRaRqQKjSh9uIdpHIol\n/0L3N20qvWziuHGwyy7hNnlyuOW3Sz02+rn29uo/j8pFcagsI1KFQsmopyfZsfdmoaW9007Q2VnZ\nvxkaCtcF3H47nHdeWF3rjjvCNQJ77BE6nDdvHv65eTP84Q9hwZb8Y4OD5d+no6P6E8QOO4RvElde\nCZ/4REjqJ58cv1ykk06g5C6SYm1t4YKw888P/QC1dEJv2TKc+AudDEY/lv/50kvD26+9VrzD+dOf\nHt4ePx4+9alwAtt55+GfxbYrfayaDmr1UQRK7iIpN9ZO6IkTw0Viu+1WewzbtoXJ4gqdFH7wg/Bt\nYs6csFbA66+H4aivvz68/fzzw4/lfw4NVf7+7e3VnTBOOgmOOw6OOALuvhsuvjgsTvPsszBpUrhN\nnNjY/opGf6NQzV1EapZvFZ9zTrhmodJvE+7hG8XohD9yu9rH3ngjfMOo9KQxfvxwoq/ltssuw9s7\n7xxer5JjNdZhwKq5i0hdjU5OPT2VJyuzUOfv6Eh2tS/30Mdw111hWOuJJ4bJ6i68MExi99pr5W8v\nvLD9/WJXSY+2ww7lTwhHHQV/+7dwyimhzFbPOZWU3EWkJo24ZqFaZvDQQ6HDecWK4RNO/qRzwgnV\nvV7+G0apk8GrrxZ/7qWXYN264fv5dRCWLg1XaNfzOKksIyJNJc2jZe69N3RAz58PV11VW8td49xF\nRFKk0TX3Ihc8i4hIkkqVsepBLXcRkQxJfVnGzAaA39b4z6cCGxMMJymKqzqKq3ppjU1xVWcscb3L\n3aeV2ylach8LM+uv5MzVaIqrOoqremmNTXFVpxFxqeYuItKElNxFRJpQVpP70tgBFKG4qqO4qpfW\n2BRXdeoeVyZr7iIiUlpWW+4iIlJCppK7mV1rZhvMbE3sWEYys73MbKWZrTWzp8xsQeyYAMysw8x+\nZmaP5+L6euyYRjKzNjP7uZndHjuWPDN7zsyeNLPHzCw1F2KY2a5mtsLMfpn7PZudgpjenztO+dur\nZnZe7LgAzOz/5H7n15jZDWbWETsmADNbkIvpqXofq0yVZczsUOB1YJm7fyR2PHlmNh2Y7u6Pmtkk\nYDVwvLv/InJcBuzk7q+b2QTgQWCBuz8cM648M/sS0AXs4u7Hxo4HQnIHutw9VWOjzew64AF3v9rM\n2oEd3f2V2HHlmVkb8DzwMXev9fqVpGLZk/C7/iF3/5OZLQfudPfvR47rI8CNwEHAIHAXcI67/6oe\n75eplru73w+8HDuO0dz9RXd/NLf9GrAW2DNuVOBBbh46JuRuqTibm9kM4Bjg6tixpJ2Z7QIcClwD\n4O6DaUrsOXOBX8dO7COMB3Yws/HAjsALkeMB+CDwsLu/6e5bgfuA/1WvN8tUcs8CM+sEZgGPxI0k\nyJU+HgM2AHe7eyriApYAC4ESyz9H4cB/m9lqM5sfO5icvwQGgO/lylhXm9lOsYMa5STghthBALj7\n88A3gd8BLwKb3f2/40YFwBrgUDObYmY7AkcDe9XrzZTcE2RmOwM3A+e5+6ux4wFw9yF33w+YARyU\n+2oYlZkdC2xw99WxYyngEHffHzgK+EKuFBjbeGB/4Ap3nwW8AVwQN6RhuTLRccD/jR0LgJm9A5gH\nvBt4J7CTmZ0aNypw97XAN4C7CSWZx4EKlwKpnpJ7QnI17ZuB6939ltjxjJb7Gt8LHBk5FIBDgONy\n9e0bgcPM7IdxQwrc/YXczw3AfxHqo7GtB9aP+Na1gpDs0+Io4FF3fyl2IDl/A/zG3Qfc/W3gFuDj\nkWMCwN2vcff93f1QQom5LvV2UHJPRK7j8hpgrbt/O3Y8eWY2zcx2zW3vQPil/2XcqMDdF7n7DHfv\nJHyd/4m7R29ZmdlOuQ5xcmWPTxC+Skfl7n8Afm9m7889NBeI2lk/ymdISUkm53fAwWa2Y+5vcy6h\nHyw6M9st93Nv4ATqeNwytcyemd0AdANTzWw9cJG7XxM3KiC0RE8DnszVtwEudPc7I8YEMB24LjeS\nYRyw3N1TM+wwhXYH/ivkA8YD/+nud8UN6X98Ebg+VwJ5Fvhs5HgAyNWODwfOjh1Lnrs/YmYrgEcJ\nZY+fk54rVW82synA28AX3P2P9XqjTA2FFBGRyqgsIyLShJTcRUSakJK7iEgTUnIXEWlCSu4iIk1I\nyV1EpAkpuYuINCEldxGRJvT/AQay1mRozexKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x297eb13b198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.grid(True)\n",
    "plt1 = plt.subplot(2, 1, 1)\n",
    "plt1.plot(x[:,0], x[:, 1], 'k.')\n",
    "plt2 = plt.subplot(2, 1, 2)\n",
    "plt2.plot(k, meandistortions, 'bx-') # 根据成本自己选个K吧，貌似4更好好些，9的话粒度太细了。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
